# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <math.h>
# include <stdlib.h>

# ifdef _OMP_VEL_ADVANCE_
# include <omp.h>
# endif

extern long int solve_poisson(struct grid *g, double ****rhs, unsigned short iRK);

extern void vel_bc(struct grid *g, short tflag, unsigned short iRK);
extern void press_bc(struct grid *g, short tflag, unsigned short iRK);

extern double dens(double T);
extern double kvisc(double T);
extern short interpol_2dim(struct var_state *var, struct cmtrx *cmx, struct cmtrx *cmy, short icellx, short icelly, short icellz);

static double getA(int i, int j, int k, struct grid *g, unsigned short isec);
static double getB(int i, int j, int k, struct grid *g, unsigned short isec);
static double getC(int i, int j, int k, struct grid *g, unsigned short isec);

struct diverg_max
     {
      short im, jm, km;
      double divmax;
     };

static void divergence(struct grid *g, double ***diverg, struct diverg_max *mdiv, unsigned short isec);

extern unsigned int vel_advance(struct grid *g, unsigned int *pcnt_max, struct arr_struct *depVar, struct arr_struct (*memFace)[3], int t_cnt)
{
 unsigned int pcnt, pcnt_subttl;

 short i, j, k, _i, _j, _k, nr, nth, nz, isec;
 double density0, dt, *r, *th, *z, ***d[Nsec], ***D0[Nsec],\
        ***_u[Nsec][nRK+1], ***_v[Nsec][nRK+1], ***_w[Nsec][nRK+1], ***_p[Nsec][nRK],\
        ***p, ***u0, ***v0, ***w0, ***p0, ***u1, ***v1, ***w1, ***p1, ***u2, ***v2, ***w2, ***p2, ***u3, ***v3, ***w3,\
        ***_A[Nsec][nRK], ***_B[Nsec][nRK], ***_C[Nsec][nRK], ***A0, ***B0, ***C0, ***A1, ***B1, ***C1, ***A2, ***B2, ***C2,\
        ***divH[Nsec][nRK-1], ***divH0, ***divH1, divH2_ijk, ***div2p0[Nsec], div2p1_ijk;

// struct timeb t1, t2;

 th = g->th; dt = g->dt; nth = g->nth;

 *pcnt_max = pcnt_subttl = 0;
 density0 = g->density0;

/* Memory for all section of the computational domain. */
  for(isec = 0; isec < Nsec; isec++)
    {
     _u[isec][0] = g->u[isec];         _v[isec][0] = g->v[isec];         _w[isec][0] = g->w[isec];          _p[isec][0] = g->p[isec];
     _u[isec][1] = depVar[isec].q[4];  _v[isec][1] = depVar[isec].q[5];  _w[isec][1] = depVar[isec].q[6];   _p[isec][1] = depVar[isec].q[7];
     _u[isec][2] = depVar[isec].q[8];  _v[isec][2] = depVar[isec].q[9];  _w[isec][2] = depVar[isec].q[10];  _p[isec][2] = depVar[isec].q[11];
     _u[isec][3] = depVar[isec].q[12]; _v[isec][3] = depVar[isec].q[13]; _w[isec][3] = depVar[isec].q[14];

     D0[isec] = depVar[isec].q[11]; divH[isec][0] = depVar[isec].q[12]; divH[isec][1] = depVar[isec].q[13]; div2p0[isec] = depVar[isec].q[14];

     for(i = 0; i < nRK; i++)
         {
          _A[isec][i] = memFace[isec][0].q[i];
          _B[isec][i] = memFace[isec][1].q[i];
          _C[isec][i] = memFace[isec][2].q[i];
         }
    }

/* Calculating divergence field at the current n-th time-step. It gets filled in D0[isec][][][]. */
  for(isec = 0; isec < Nsec; isec++) divergence(g, D0[isec], (struct diverg_max *)0, isec);

/* --------------------------- Calculation of u(1), v(1), w(1), and p(1) ---------------------------------- */
  for(isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0];
       u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0]; p0 = _p[isec][0];
       u1 = _u[isec][1]; v1 = _v[isec][1]; w1 = _w[isec][1];


#      ifdef _OMP_VEL_ADVANCE_
       #pragma omp parallel private(i,j,k)
#      endif
          {
#          ifdef _OMP_VEL_ADVANCE_
            #pragma omp for schedule(dynamic,1)
#          endif
           for(k = 0; k <= nz; k++) for(i = 0; i <= nr; i++) for(j = 0; j <= nth; j++)
              {
               if(j < nth && k < nz) A0[i][j][k] = getA(NrG+i-1, NthG+j, NzG+k, g, isec);

               if(i < nr && k < nz) B0[i][j][k] = getB(NrG+i, NthG+j-1, NzG+k, g, isec);
               if(i < nr && j < nth) C0[i][j][k] = getC(NrG+i, NthG+j, NzG+k-1, g, isec);
              }

#          ifdef _OMP_VEL_ADVANCE_
            #pragma omp for schedule(dynamic,1)
#          endif
           for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
              {
               u1[i][j][k] = u0[i][j][k]\
                 + c2*dt*( A0[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p0[i+1][j][k] - p0[i][j][k])/(r[i+2]-r[i]) );

               v1[i][j][k] = v0[i][j][k]\
                 + c2*dt*( B0[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p0[i][j+1][k] - p0[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) );

               w1[i][j][k] = w0[i][j][k]\
                 + c2*dt*( C0[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p0[i][j][k+1] - p0[i][j][k])/(z[k+2]-z[k]) );

               if(i == NrG+nr-1) u1[i][j][k] = 0;
               if(isec == 1 && k == NzG+nz-1) w1[i][j][k] = 0;
              }
          }

       g->u[isec] = u1; g->v[isec] = v1; g->w[isec] = w1;
      }

  vel_bc(g, 0, 1);

/* _u[1], _v[1], _w[1] is used in "press_bc()". Therefore, the memory used by these can't be distroyed. */
  for(isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0]; p0 = _p[isec][0]; divH0 = divH[isec][0];
       A1 = _A[isec][1]; B1 = _B[isec][1]; C1 = _C[isec][1]; p1 = _p[isec][1]; divH1 = divH[isec][1];

       d[isec] = depVar[isec].q[8]; p = _p[isec][0];

#      ifdef _OMP_VEL_ADVANCE_
       #pragma omp parallel private(i,j,k,_i,_j,_k)
#      endif
         {
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic,1)
#         endif
          for(k = 0; k <= nz; k++) for(i = 0; i <= nr; i++) for(j = 0; j <= nth; j++)
              {
               if(j < nth && k < nz) A1[i][j][k] = getA(NrG+i-1, NthG+j, NzG+k, g, isec);
               if(i < nr && k < nz)  B1[i][j][k] = getB(NrG+i, NthG+j-1, NzG+k, g, isec);
               if(i < nr && j < nth) C1[i][j][k] = getC(NrG+i, NthG+j, NzG+k-1, g, isec);
              }

       /* Calculating RHS of the PPE-1. */
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic,1) nowait
#         endif
          for(_k = 0; _k < nz; _k++) for(k = _k+NzG, _i = 0; _i < nr; _i++) for(i = _i+NrG, _j = 0, j = NthG; _j < nth; j++, _j++)
               {
                div2p0[isec][i][j][k] = DIV2P_ijk;

                divH0[i][j][k] = 2*(( r[i+1]*A0[_i+1][_j][_k] - ((((int)(r[i]*10000)) == 0) ? A0[_i][_j][_k]:r[i]*A0[_i][_j][_k]))/(r[i+1]-r[i])\
                           + (B0[_i][_j+1][_k] - B0[_i][_j][_k])/(th[j+1]-th[j]) )/(r[i]+r[i+1])\
                           + (C0[_i][_j][_k+1] - C0[_i][_j][_k])/(z[k+1]-z[k]);

                divH1[i][j][k] = 2*( (r[i+1]*A1[_i+1][_j][_k] - ((((int)(r[i]*10000)) == 0) ? A1[_i][_j][_k]:r[i]*A1[_i][_j][_k]))/(r[i+1]-r[i])\
                          + (B1[_i][_j+1][_k] - B1[_i][_j][_k])/(th[j+1]-th[j]) )/(r[i]+r[i+1])\
                          + (C1[_i][_j][_k+1] - C1[_i][_j][_k])/(z[k+1]-z[k]);

                d[isec][i][j][k] = density0*( divH1[i][j][k]\
                                            + (D0[isec][i][j][k]/dt + a31*(divH0[i][j][k] - div2p0[isec][i][j][k]/density0))/a32 );
               }

       /* Copying _p[0] into _p[1] to provide it as initial condition for pressure while solving the PPE. */
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic,2)
#         endif
          for(k = 0; k < 2*NzG+nz; k++) for(i = 0; i < 2*NrG+nr; i++) for(j = 0; j < 2*NthG+nth; j++)
             p1[i][j][k] = p0[i][j][k];
         }
       g->p[isec] = p1; /* Used in solve_poisson(). */
       g->A[isec] = A1; /* Used in press_bc(). */
       g->B[isec] = B1; /* Used in press_bc(). */
       g->C[isec] = C1; /* Used in press_bc(). */
      }

/* Solving the PPE(1) to calculate p(1). */
   if(g->tptr != (FILE *)0) fprintf(g->tptr, "(PPE - 1)\n");
   pcnt = solve_poisson(g, d, 1);

   pcnt_subttl += pcnt;
   if(pcnt > *pcnt_max) *pcnt_max = pcnt;
/* ------------------------------------------------------------------------------------------------------------------------- */


/* ------------------------------------- Calculation of u(2), v(2), w(2), and p(2) ----------------------------------------- */
/* Calculating u(2), v(2), and w(2). */

  for(isec = 0; isec < Nsec; isec++)
     {
      nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec]; p0 = _p[isec][0];  p1 = _p[isec][1];
      u2 = _u[isec][2]; v2 = _v[isec][2]; w2 = _w[isec][2]; u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0];
      A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0]; A1 = _A[isec][1]; B1 = _B[isec][1]; C1 = _C[isec][1];

#     ifdef _OMP_VEL_ADVANCE_
      #pragma omp parallel for private(i,j,k) schedule(dynamic,2)
#     endif
      for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
          {
           u2[i][j][k] = u0[i][j][k]\
                + a31*dt*( A0[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p0[i+1][j][k] - p0[i][j][k])/(r[i+2] - r[i]) )\
                + a32*dt*( A1[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p1[i+1][j][k] - p1[i][j][k])/(r[i+2] - r[i]) );

           v2[i][j][k] = v0[i][j][k]\
                + a31*dt*( B0[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p0[i][j+1][k] - p0[i][j][k])/( (r[i]+r[i+1])*(th[j+2]-th[j]) ) )\
                + a32*dt*( B1[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p1[i][j+1][k] - p1[i][j][k])/( (r[i]+r[i+1])*(th[j+2]-th[j]) ) );

           w2[i][j][k] = w0[i][j][k]\
                + a31*dt*( C0[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p0[i][j][k+1] - p0[i][j][k])/(z[k+2]-z[k]) )\
                + a32*dt*( C1[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p1[i][j][k+1] - p1[i][j][k])/(z[k+2]-z[k]) );

           if(i == NrG+nr-1) u2[i][j][k] = 0;
           if(isec == 1 && k == NzG+nz-1) w2[i][j][k] = 0;
          }

      g->u[isec] = u2; g->v[isec] = v2; g->w[isec] = w2;
     }


  vel_bc(g, 0, 2);


/* _u[2], _v[2], _w[2] is used in "press_bc()". Therefore, the memory used by these can't be distroyed. */

  for(isec = 0; isec < Nsec; isec++)
     {
      nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
      A2 = _A[isec][2]; B2 = _B[isec][2]; C2 = _C[isec][2]; p1 = _p[isec][1]; p2 = _p[isec][2];
      d[isec] = depVar[isec].q[4]; p = _p[isec][1]; divH0 = divH[isec][0]; divH1 = divH[isec][1];

#     ifdef _OMP_VEL_ADVANCE_
      #pragma omp parallel private(i,j,k,_i,_j,_k,div2p1_ijk,divH2_ijk)
#     endif
          {
#          ifdef _OMP_VEL_ADVANCE_
           #pragma omp for schedule(dynamic,1)
#          endif
           for(k = 0; k <= nz; k++) for(i = 0; i <= nr; i++) for(j = 0; j <= nth; j++)
               {
                if(j < nth && k < nz) A2[i][j][k] = getA(NrG+i-1, NthG+j, NzG+k, g, isec);
                if(i < nr && k < nz)  B2[i][j][k] = getB(NrG+i, NthG+j-1, NzG+k, g, isec);
                if(i < nr && j < nth) C2[i][j][k] = getC(NrG+i, NthG+j, NzG+k-1, g, isec);
               }

#          ifdef _OMP_VEL_ADVANCE_
           #pragma omp for schedule(dynamic,1)
#          endif
           for(_k = 0; _k < nz; _k++) for(k = _k+NzG, _i = 0; _i < nr; _i++) for(i = _i+NrG, _j = 0, j = NthG; _j < nth; j++, _j++)
               {
                div2p1_ijk = DIV2P_ijk;
                divH2_ijk = 2*( (r[i+1]*A2[_i+1][_j][_k] - ( (((int)(r[i]*10000)) == 0) ? A2[_i][_j][_k]:r[i]*A2[_i][_j][_k]))/(r[i+1]-r[i])\
                     + (B2[_i][_j+1][_k] - B2[_i][_j][_k])/(th[j+1]-th[j]) )/(r[i]+r[i+1])\
                     + (C2[_i][_j][_k+1] - C2[_i][_j][_k])/(z[k+1]-z[k]);

                d[isec][i][j][k] = density0*( divH2_ijk + (D0[isec][i][j][k]/dt + a41*(divH0[i][j][k] - div2p0[isec][i][j][k]/density0)\
                                                         + a42*(divH1[i][j][k] - div2p1_ijk/density0))/a43 );
               }

        /* Copying _p[1] into _p[2] to provide it as initial condition for pressure while solving the PPE. */
#          ifdef _OMP_VEL_ADVANCE_
           #pragma omp for schedule(static)
#          endif
        /* IMPORTANT: The "D0" and "p2" are sharing same memory space i.e. the above calculation of "d" must be over before "p2"
           gets written below. i.e. we can't use the "nowait" flag in the above parallelization. */
           for(k = 0; k < 2*NzG+nz; k++) for(i = 0; i < 2*NrG+nr; i++) for(j = 0; j < 2*NthG+nth; j++)
            p2[i][j][k] = p1[i][j][k];
          }
      g->p[isec] = p2;
      g->A[isec] = A2; /* Used in press_bc(). */
      g->B[isec] = B2; /* Used in press_bc(). */
      g->C[isec] = C2; /* Used in press_bc(). */
     }


/* Solving the PPE(2) to calculate p(2). */
   if(g->tptr != (FILE *)0) fprintf(g->tptr, "(PPE - 2)\n");
   pcnt = solve_poisson(g, d, 2);

   pcnt_subttl += pcnt;
   if(pcnt > *pcnt_max) *pcnt_max = pcnt;
/* ---------------------------------------------------------------------------------------------------------- */

/* ================================= Calculation of u, v, w, and p at next time-step ========================== */
/* Calculation of u, v, and w at next time-step. */
  for(isec = 0; isec < Nsec; isec++)
      {
       nr = g->nr[isec]; nz = g->nz[isec]; r = g->r[isec]; z = g->z[isec];
       u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0];
       u3 = _u[isec][3]; v3 = _v[isec][3]; w3 = _w[isec][3];
       A0 = _A[isec][0]; B0 = _B[isec][0]; C0 = _C[isec][0]; p0 = _p[isec][0];
       A1 = _A[isec][1]; B1 = _B[isec][1]; C1 = _C[isec][1]; p1 = _p[isec][1];
       A2 = _A[isec][2]; B2 = _B[isec][2]; C2 = _C[isec][2]; p2 = _p[isec][2];

#      ifdef _OMP_VEL_ADVANCE_
       #pragma omp parallel for private(i,j,k) schedule(dynamic,2)
#      endif
       for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
           {
            u3[i][j][k] = u0[i][j][k]\
                  + b1*dt*( A0[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p0[i+1][j][k] - p0[i][j][k])/(r[i+2]-r[i]) )\
                  + b2*dt*( A1[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p1[i+1][j][k] - p1[i][j][k])/(r[i+2]-r[i]) )\
                  + b3*dt*( A2[i-NrG+1][j-NthG][k-NzG] - (2.0/density0)*(p2[i+1][j][k] - p2[i][j][k])/(r[i+2]-r[i]) );

            v3[i][j][k] = v0[i][j][k]\
                  + b1*dt*( B0[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p0[i][j+1][k] - p0[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) )\
                  + b2*dt*( B1[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p1[i][j+1][k] - p1[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) )\
                  + b3*dt*( B2[i-NrG][j-NthG+1][k-NzG] - (4.0/density0)*(p2[i][j+1][k] - p2[i][j][k])/((r[i]+r[i+1])*(th[j+2]-th[j])) );

            w3[i][j][k] = w0[i][j][k]\
                  + b1*dt*( C0[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p0[i][j][k+1] - p0[i][j][k])/(z[k+2]-z[k]) )\
                  + b2*dt*( C1[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p1[i][j][k+1] - p1[i][j][k])/(z[k+2]-z[k]) )\
                  + b3*dt*( C2[i-NrG][j-NthG][k-NzG+1] - (2.0/density0)*(p2[i][j][k+1] - p2[i][j][k])/(z[k+2]-z[k]) );

            if(i == NrG+nr-1) u3[i][j][k] = 0;
            if(isec == 1 && k == NzG+nz-1) w3[i][j][k] = 0;
           }

       g->u[isec] = u3; g->v[isec] = v3; g->w[isec] = w3;
      }

  vel_bc(g, 1, 0);

/* Storing mean velocities and copying the newly evaluated velocities to _u[0], _v[0], and _w[0]. */
  for(isec = 0; isec < Nsec; isec++)
     {
      nr = g->nr[isec]; nz = g->nz[isec];
      p0 = _p[isec][0]; p1 = _p[isec][1]; p2 = _p[isec][2];
      u0 = _u[isec][0]; v0 = _v[isec][0]; w0 = _w[isec][0];
      u3 = _u[isec][3]; v3 = _v[isec][3]; w3 = _w[isec][3];

   /* Memory to store mean velocities at time = t + dt/2. */
      g->um[isec] = depVar[isec].q[4]; g->vm[isec] = depVar[isec].q[5]; g->wm[isec] = depVar[isec].q[6];

#     ifdef _OMP_VEL_ADVANCE_
      #pragma omp parallel private(i,j,k)
#     endif
         {
#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(static) nowait
#         endif
          for(k = 0; k < nz+2*NzG; k++) for(i = 0; i < nr+2*NrG; i++) for(j = 0; j < nth+2*NthG; j++)
              {
            /* Storing mean-velocities. */
               g->um[isec][i][j][k] = 0.5*(u0[i][j][k] + u3[i][j][k]);
               g->vm[isec][i][j][k] = 0.5*(v0[i][j][k] + v3[i][j][k]);
               g->wm[isec][i][j][k] = 0.5*(w0[i][j][k] + w3[i][j][k]);

            /* Copying newly evaluated velocities. */
               u0[i][j][k] = u3[i][j][k];
               v0[i][j][k] = v3[i][j][k];
               w0[i][j][k] = w3[i][j][k];
              }

#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(static)
#         endif
          for(k = NzG; k < nz + NzG; k++) for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
          /* Pressure at the next time-step will be given by _p[0]. */
             p0[i][j][k] = d1*p0[i][j][k] + d2*p1[i][j][k] + d3*p2[i][j][k];
         }
      g->u[isec] = u0; g->v[isec] = v0; g->w[isec] = w0; g->p[isec] = p0;
     }

  press_bc(g, 1, 0);
/* ---------------------------------------------------------------------------------------------------------------- */

 return pcnt_subttl;

/*** q[0]~q[6] are engagged. These can't be used in other functions. ***/
}


/* For (i,j,k)-th CELL of isec-th section. */
static double getA(int i, int j, int k, struct grid *g, unsigned short isec)
  {
   double ***u, ***v, ***w, ***T, *r, *th, *z, kvisc0;
   struct cmtrx *cm, *cm2;
   double q, visc_term, conv_term;
   unsigned short indx;
   struct var_state var;

   u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; T = g->T[isec]; r = g->r[isec];  th = g->th;  z = g->z[isec];
   kvisc0 = g->kvisc0;

#  ifndef DISCR_AUTO
   return ((int)(r[i+1]*10000) == 0) ? rA_CENTER : AA;
#  endif

   visc_term = conv_term = 0.0;

/* Calculating u*(du/dr) at r-face centers. */
// Evaluating du/dr at r-face center.
   cm = &g->dcmx_rss1[isec];
   for(q = (double)0, indx = 0; indx <= cm->iend[i+1-cm->respadi]-cm->istart[i+1-cm->respadi]; indx++)
   q += cm->mtrx[i+1-cm->respadi][indx]*u[cm->srcpadi-1+cm->istart[i+1-cm->respadi]+indx][j][k];

/* If r[i+1] == 0 (where i = cell no.), return [r*du/dr]_FaceCenter because other terms become zero. */
   if(((int)(r[i+1]*10000)) == 0) return kvisc0*q;

   conv_term += q*u[i][j][k];

/* Calculating (1/r)*(du/dr) : a viscous term. */
   visc_term += q/r[i+1];

/* Calculating (v/r)*(du/dth) at r-face centers. */
// Determining "v" at r-face center.
   cm = &g->icmx_rcs[isec];
   cm2 = &g->icmx_thsc;
   var.q = v;
   var.dirx = 1;
   var.diry = 2;
   var.dxdef = 0.0;
   var.dydef = 0.5;
   var.dxreq = 0.5;
   var.dyreq = 0.0;
   interpol_2dim(&var, cm, cm2, i, j, k);

// Calculating du/dth at r-face center.
  cm = &g->dcmx_thcc1;
  for(q = 0.0, indx = 0; indx <= cm->iend[j-cm->respadi]-cm->istart[j-cm->respadi]; indx++)
   q += cm->mtrx[j-cm->respadi][indx]*u[i][cm->srcpadi+cm->istart[j-cm->respadi]+indx][k];

  conv_term += var.qinterp*q/r[i+1];

/* Calculating v^2/r at the r-face center. */
  conv_term -= var.qinterp*var.qinterp/r[i+1];

/* Calculating w*du/dz at r-face center. */
// Determining "w" at r-face center.
  cm = &g->icmx_rcs[isec];
  cm2 = &g->icmx_zsc[isec][0];
  if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm2 = &g->icmx_zsc[isec][1];
  var.q = w;
  var.dirx = 1;
  var.diry = 3;
  var.dxdef = 0.0;
  var.dydef = 0.5;
  var.dxreq = 0.5;
  var.dyreq = 0.0;
  interpol_2dim(&var, cm, cm2, i, k, j);

// Calculating du/dz at r-face center.
  cm = &g->dcmx_zcc1[isec][0];
  if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zcc1[isec][1];
  for(q = 0.0, indx = 0; indx <= cm->iend[k-cm->respadi]-cm->istart[k-cm->respadi]; indx++)
   q += cm->mtrx[k-cm->respadi][indx]*u[i][j][cm->srcpadi+cm->istart[k-cm->respadi]+indx];

  conv_term += var.qinterp*q;

/* Calculating d2u/dr2 at r-face centers: a viscous term */
  cm = &g->dcmx_rss2[isec];
  for(q = (double)0, indx = 0; indx <= cm->iend[i+1-cm->respadi]-cm->istart[i+1-cm->respadi]; indx++)
  q += cm->mtrx[i+1-cm->respadi][indx]*u[cm->srcpadi-1+cm->istart[i+1-cm->respadi]+indx][j][k];

  visc_term += q;

/* Calculating (1/r^2)*(d2u/dth2) at r-face centers: a viscous term */
  cm = &g->dcmx_thcc2;
  for(q = 0.0, indx = 0; indx <= cm->iend[j-cm->respadi]-cm->istart[j-cm->respadi]; indx++)
   q += cm->mtrx[j-cm->respadi][indx]*u[i][cm->srcpadi+cm->istart[j-cm->respadi]+indx][k];

  visc_term += (q/r[i+1])/r[i+1];

/* Calculating d2u/dz2 at r-face center: a viscous term. */
  cm = &g->dcmx_zcc2[isec][0];
  if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zcc2[isec][1];
  for(q = 0.0, indx = 0; indx <= cm->iend[k-cm->respadi]-cm->istart[k-cm->respadi]; indx++)
   q += cm->mtrx[k-cm->respadi][indx]*u[i][j][cm->srcpadi+cm->istart[k-cm->respadi]+indx];

  visc_term += q;

/* Calculating (1/r^2)*(u + 2*dv/dth) at r-face center: a viscous term. */
// Calculating dv/dth at r-face center.
  cm = &g->icmx_rcs[isec];
  cm2 = &g->dcmx_thsc1;
  var.q = v;
  var.dirx = 1;
  var.diry = 2;
  var.dxdef = 0.0;
  var.dydef = 0.5;
  var.dxreq = 0.5;
  var.dyreq = 0.0;
  interpol_2dim(&var, cm, cm2, i, j, k);

  visc_term -= ((u[i][j][k]+2*var.qinterp)/r[i+1])/r[i+1];

  return kvisc0*visc_term - conv_term;
 }


/* For (i,j,k)-th CELL of isec-th section. */
static double getB(int i, int j, int k, struct grid *g, unsigned short isec)
  {
   double ***u, ***v, ***w, ***T, *r, *th, *z, kvisc0;
   struct cmtrx *cm, *cm2;
   double q, visc_term, conv_term;
   unsigned short indx;
   struct var_state var;

   u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; T = g->T[isec];  r = g->r[isec];  th = g->th;  z = g->z[isec];
   kvisc0 = g->kvisc0;

#  ifndef DISCR_AUTO
   return BB;
#  endif

   conv_term = visc_term = 0.0;

/* Evaluating u.dv/dr at the th-face center. */
// Calculating u at th-face center.
   cm = &g->icmx_rsc[isec];
   cm2 = &g->icmx_thcs;
   var.q = u;
   var.dirx = 1;
   var.diry = 2;
   var.dxdef = 0.5;
   var.dydef = 0.0;
   var.dxreq = 0.0;
   var.dyreq = 0.5;
   interpol_2dim(&var, cm, cm2, i, j, k);

/* Evaluating u*v/r at th-face center. */
   conv_term += var.qinterp*2*v[i][j][k]/(r[i]+r[i+1]);

// Calculating dv/dr at th-face center.
   cm = &g->dcmx_rcc1[isec];
   for(q = 0.0, indx = 0; indx <= cm->iend[i-cm->respadi]-cm->istart[i-cm->respadi]; indx++)
      q += cm->mtrx[i-cm->respadi][indx]*v[cm->srcpadi+cm->istart[i-cm->respadi]+indx][j][k];

   conv_term += var.qinterp*q;

/* Evaluating (1/r)*(dv/dr) at th-face center. */
   visc_term += 2*q/(r[i]+r[i+1]);

/* Evaluating (v/r)*(dv/dth) at th-face center. */
// Calculating dv/dth at th-face center.
   cm = &g->dcmx_thss1;
   for(q = 0.0, indx = 0; indx <= cm->iend[j+1-cm->respadi]-cm->istart[j+1-cm->respadi]; indx++)
     q += cm->mtrx[j+1-cm->respadi][indx]*v[i][cm->srcpadi-1+cm->istart[j+1-cm->respadi]+indx][k];

   conv_term += 2*q*v[i][j][k]/(r[i]+r[i+1]);

/* Evaluating w*dv/dz at the th-face center. */
// Calculating w at th-face center.
   cm = &g->icmx_thcs;
   cm2 = &g->icmx_zsc[isec][0];
   if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm2 = &g->icmx_zsc[isec][1];
   var.q = w;
   var.dirx = 2;
   var.diry = 3;
   var.dxdef = 0.0;
   var.dydef = 0.5;
   var.dxreq = 0.5;
   var.dyreq = 0.0;
   interpol_2dim(&var, cm, cm2, j, k, i);

// Calculating dv/dz at th-face center.
   cm = &g->dcmx_zcc1[isec][0];
   if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zcc1[isec][1];
   for(q = 0.0, indx = 0; indx <= cm->iend[k-cm->respadi]-cm->istart[k-cm->respadi]; indx++)
     q += cm->mtrx[k-cm->respadi][indx]*v[i][j][cm->srcpadi+cm->istart[k-cm->respadi]+indx];

   conv_term += var.qinterp*q;

/* Evaluating d2v/dr2 at the th-face center. */
   cm = &g->dcmx_rcc2[isec];
   for(q = 0.0, indx = 0; indx <= cm->iend[i-cm->respadi]-cm->istart[i-cm->respadi]; indx++)
      q += cm->mtrx[i-cm->respadi][indx]*v[cm->srcpadi+cm->istart[i-cm->respadi]+indx][j][k];

   visc_term += q;

/* Evaluating (1/r^2)*(d2v/dth2) at th-face center. */
// Calculating d2v/dth2 at th-face center.
   cm = &g->dcmx_thss2;
   for(q = 0.0, indx = 0; indx <= cm->iend[j+1-cm->respadi]-cm->istart[j+1-cm->respadi]; indx++)
     q += cm->mtrx[j+1-cm->respadi][indx]*v[i][cm->srcpadi-1+cm->istart[j+1-cm->respadi]+indx][k];

   visc_term += 4*(q/(r[i]+r[i+1]))/(r[i]+r[i+1]);

/* Evaluating d2v/dz2 at th-face center. */
// Calculating d2v/dz2 at th-face center.
   cm = &g->dcmx_zcc2[isec][0];
   if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zcc2[isec][1];
   for(q = 0.0, indx = 0; indx <= cm->iend[k-cm->respadi]-cm->istart[k-cm->respadi]; indx++)
     q += cm->mtrx[k-cm->respadi][indx]*v[i][j][cm->srcpadi+cm->istart[k-cm->respadi]+indx];

   visc_term += q;

/* Evaluating (v - 2*du/dth)/r^2 at th-face center. */
// Calculating du/dth at the th-face center.
   cm = &g->icmx_rsc[isec];
   cm2 = &g->dcmx_thcs1;
   var.q = u;
   var.dirx = 1;
   var.diry = 2;
   var.dxdef = 0.5;
   var.dydef = 0.0;
   var.dxreq = 0.0;
   var.dyreq = 0.5;
   interpol_2dim(&var, cm, cm2, i, j, k);

   visc_term -= 4*((v[i][j][k] - 2*var.qinterp)/(r[i]+r[i+1]))/(r[i]+r[i+1]);

   return kvisc0*visc_term - conv_term;
  }


/* For (i,j,k)-th CELL of isec-th section. */
static double getC(int i, int j, int k, struct grid *g, unsigned short isec)
  {
   double ***u, ***v, ***w, ***T, *r, *th, *z, kvisc0, density0;
   struct cmtrx *cm, *cm2;
   double q, visc_term, conv_term;
   unsigned short indx;
   struct var_state var;

   u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; T = g->T[isec];  r = g->r[isec];  th = g->th;  z = g->z[isec];
   kvisc0 = g->kvisc0; density0 = g->density0;

#  ifndef DISCR_AUTO
   return CC;
#  endif

   visc_term = conv_term = 0.0;

/* Evaluating u*dw/dr at z-face center. */
// Calculating u at z-face center.
   cm = &g->icmx_rsc[isec];
   cm2 = &g->icmx_zcs[isec][0];
   if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm2 = &g->icmx_zcs[isec][1];
   var.q = u;
   var.dirx = 1;
   var.diry = 3;
   var.dxdef = 0.5;
   var.dydef = 0.0;
   var.dxreq = 0.0;
   var.dyreq = 0.5;
   interpol_2dim(&var, cm, cm2, i, k, j);

// Calculating dw/dr at z-face center.
   cm = &g->dcmx_rcc1[isec];
   for(q = 0.0, indx = 0; indx <= cm->iend[i-cm->respadi]-cm->istart[i-cm->respadi]; indx++)
     q += cm->mtrx[i-cm->respadi][indx]*w[cm->srcpadi+cm->istart[i-cm->respadi]+indx][j][k];

   conv_term += var.qinterp*q;

/* Evaluating (1/r)*(dw/dr) at w-face center. */
   visc_term += 2*q/(r[i]+r[i+1]);


/* Evaluating (v/r)*(dw/dth) at z-face center. */

// Calculating v at w-face center.
   cm = &g->icmx_thsc;
   cm2 = &g->icmx_zcs[isec];
   var.q = v;
   var.dirx = 2;
   var.diry = 3;
   var.dxdef = 0.5;
   var.dydef = 0.0;
   var.dxreq = 0.0;
   var.dyreq = 0.5;
   interpol_2dim(&var, cm, cm2, j, k, i);

// Calculating (dw/dth) at z-face center.
   cm = &g->dcmx_thcc1;
   for(q = 0.0, indx = 0; indx <= cm->iend[j-cm->respadi]-cm->istart[j-cm->respadi]; indx++)
     q += cm->mtrx[j-cm->respadi][indx]*w[i][cm->srcpadi+cm->istart[j-cm->respadi]+indx][k];

   conv_term += 2*var.qinterp*q/(r[i]+r[i+1]);

/* Evaluating w*dw/dz at w-face center. */
// Calculating (dw/dz) at z-face center.
   cm = &g->dcmx_zss1[isec][0];
   if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zss1[isec][1];
   for(q = 0.0, indx = 0; indx <= cm->iend[k+1-cm->respadi]-cm->istart[k+1-cm->respadi]; indx++)
     q += cm->mtrx[k+1-cm->respadi][indx]*w[i][j][cm->srcpadi-1+cm->istart[k+1-cm->respadi]+indx];

   conv_term += q*w[i][j][k];

/* Evaluating d2w/dr2 at w-face center. */
// Calculating d2w/dr2 at z-face center.
   cm = &g->dcmx_rcc2[isec];
   for(q = 0.0, indx = 0; indx <= cm->iend[i-cm->respadi]-cm->istart[i-cm->respadi]; indx++)
     q += cm->mtrx[i-cm->respadi][indx]*w[cm->srcpadi+cm->istart[i-cm->respadi]+indx][j][k];

   visc_term += q;

/* Evaluating (1/r^2)*(d2w/dth2) at z-face center. */
// Calculating d2w/dth2 at z-face center.
   cm = &g->dcmx_thcc2;
   for(q = 0.0, indx = 0; indx <= cm->iend[j-cm->respadi]-cm->istart[j-cm->respadi]; indx++)
     q += cm->mtrx[j-cm->respadi][indx]*w[i][cm->srcpadi+cm->istart[j-cm->respadi]+indx][k];

   visc_term += 4*(q/(r[i]+r[i+1]))/(r[i]+r[i+1]);

/* Evaluating d2w/dz2 at z-face center. */
// Calculating d2w/dz2 at z-face center.
   cm = &g->dcmx_zss2[isec][0];
   if(isec == 1) if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zss2[isec][1];
   for(q = 0.0, indx = 0; indx <= cm->iend[k+1-cm->respadi]-cm->istart[k+1-cm->respadi]; indx++)
     q += cm->mtrx[k+1-cm->respadi][indx]*w[i][j][cm->srcpadi-1+cm->istart[k+1-cm->respadi]+indx];

   visc_term += q;

//     printf("%0.15E\t%0.15E\n", kvisc0*visc_term - conv_term - GRAVITY_DOWNWARD, CC);

   return kvisc0*visc_term - conv_term - GRAVITY_DOWNWARD;
  }


static void divergence(struct grid *g, double ***diverg, struct diverg_max *mdiv, unsigned short isec)
   {
    double *r, *th, *z, ***u, ***v, ***w, *div;
    short i, j, k, tid, nr, nz, nth, nTRD;
    struct diverg_max *tddiv;

    u = g->u[isec];  v = g->v[isec];  w = g->w[isec]; r = g->r[isec]; th = g->th; z = g->z[isec];
    nr = g->nr[isec]; nth = g->nth; nz = g->nz[isec];

    nTRD = 1;
#   ifdef _OMP_VEL_ADVANCE_
    nTRD = g->nTRD;
#   endif

    tddiv = (struct diverg_max *)malloc(sizeof(struct diverg_max)*nTRD);
    div = (double *)malloc(sizeof(double)*nTRD);

    tid = 0;
#   ifdef _OMP_VEL_ADVANCE_
    #pragma omp parallel private(i,j,k,tid)
#   endif
         {
#         ifdef _OMP_VEL_ADVANCE_
          tid = omp_get_thread_num();
#         endif

          tddiv[tid].im = -1;
          tddiv[tid].jm = -1;
          tddiv[tid].km = -1;
          tddiv[tid].divmax = 0.0;

#         ifdef _OMP_VEL_ADVANCE_
          #pragma omp for schedule(dynamic, 3)
#         endif
          for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++) for(k = NzG; k < nz+NzG; k++)
             {
              div[tid] = DIV_ijk;

           /* Filtering out the maximum divergence. */
              if((div[tid] < 0 ? -div[tid]:div[tid]) > (tddiv[tid].divmax < 0 ? -tddiv[tid].divmax:tddiv[tid].divmax))
                 {
                  tddiv[tid].divmax = div[tid];
                  tddiv[tid].im = i;
                  tddiv[tid].jm = j;
                  tddiv[tid].km = k;
                 }

              diverg[i][j][k] = div[tid];
             }
         }

    free(div);

 /* Calculating the absolute maximum value of the divergence. */
    for(i = 0, tid = 1; tid < nTRD; tid++)
    if((tddiv[tid].divmax < 0 ? -tddiv[tid].divmax:tddiv[tid].divmax) > (tddiv[i].divmax < 0 ? -tddiv[i].divmax:tddiv[i].divmax)) i = tid;

    if(g->tptr != (FILE *)0) fprintf(g->tptr, "\tdivmax(%d, %d, %d | %d) = %E\n", tddiv[i].im, tddiv[i].jm, tddiv[i].km, isec, tddiv[i].divmax);

    if( mdiv != (struct diverg_max *)0 )
      { mdiv->divmax = tddiv[i].divmax; mdiv->im = tddiv[i].im; mdiv->jm = tddiv[i].jm; mdiv->km = tddiv[i].km; }

    free(tddiv);

    return;
   }
